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Abstract 

We demonstrate, by numerical simulations, that the dynamics of nuclear mat- 
ter mean field inside the spinodal region is chaotic. Spontaneous symmetry- 
breaking - no explicit fluctuating term is considered - occurs leading to wild 
unpredictable density fluctuations. A proper recipe to calculate an average 
Lyapunov exponent in this multidimensional phase space is introduced. The 
latter is calculated for different values of the density in order to characterize in 
a quantitative way the chaotic and regular regions. It is argued that the mean 
field chaoticity can be the main mechanism of the nuclear multifragmentation 

occurring in the intermediate energy reactions. 
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Heavy ion collisions in the intermediate energy regime between 20 and 100 MeV/nucleon 
can result in multifragment emissions, where a number of intermediate mass fragments 
(IMF : Z > 2) is observed. This disassembly of the nuclear system, which can include a 
substantial fraction of the whole colliding system, falls under the generic name of "nuclear 
multifragmentation" , and is presently under intense study in different laboratories. The 
main mechanism which is responsible of this phenomenon has not yet been determined. 
Both statistical models and dynamical simulations [|l^] have been used to describe the cluster 
formation in the multifragmentation regime. The role of mechanical instabilities, either of 
the Raleigh type 0, or associated with ring formation [|3|, has been also stressed recently. 
These calculations provide suggestions on the possible dominant path followed by the nuclear 
system along the multifragmentation process. The onset of instability in the spinodal zone 
of the nuclear matter Equation of State (EOS) has been studied in ref. j4|, by analyzing the 
corresponding linear response for imaginary frequency as a function of the wave number and 
of the density. 

In this letter we establish, by numerical simulations, that the dynamical evolution of 
nuclear matter mean field in the spinodal region is chaotic. We estimate the corresponding 
Lyapunov exponent which characterizes the rate of divergence between the trajectories in 
density space. This result can have far reaching consequences for the whole picture of the 
multifragmentation process. In particular, a) the relevance of the mean trajectory, which 
is an essential ingredient in most of the simulations with kinetic equations, appears ques- 
tionable if the spinodal region is reached; b) the onset of chaoticity implies, as we explicitly 
verify numerically, strong non-linearity of the evolution, and therefore strong coupling be- 
tween different modes; c) the occurrence of chaotic motion suggests that thermodynamical 
quantities, like entropy, are not necessarily appropriate for characterizing the system dy- 
namics, but they should be substituted by others, based on chaotic evolution concepts, like 
the Kolmogorov entropy || . Other more quantitative consequences will be discussed in the 
following. The numerical analysis is performed by solving the Vlasov equation 
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~^ + {H[f]J} = (1) 

on a two-dimensional lattice, according to the method of ref. ||. In eq. (1) H[f] is the 
self-consistent effective one-body Hamiltonian. In this approach the one-body density dis- 
tribution is represented on a grid at variance with the more commonly used test particle 
methods [0|§. The former has been preferred since it allows a better control of the numer- 
ical accuracy and conservation laws. The single particle phase space is divided into cells 
and nuclear matter is confined in a large torus with periodic boundary conditions. The side 
lengths of the torus are equal to L x = 51 fm and L y = 15 fm. In order to achieve sufficient 
accuracy, very small cells having 5x = | fm and 5p = 40 MeV/c are employed. The time 
integration step was chosen equal to 0.5 fm/c. We adopt a simplified Skyrme interaction 
for the effective one-body field which is averaged over the y-direction 0. We employ a 
Fermi momentum at saturation of Pp = 260 MeV/c, which gives a density in two dimen- 
sions equal to po = 0.55 fm~ 2 . The mean field evolution is treated by means of a standard 
matrix technique 0. The main indication of chaotic behaviour of a dynamical system is 
the extreme sensitivity of its evolution to the initial conditions. This is illustrated in figs, 
la) and lb). At the initial time nuclear matter presents a density wave, with an amplitude 
Api = 10~ 2 po ; and with a wave number k = 27r(5/L x ). The local momentum distribution 
is assumed to be the one of a Fermi gas at a temperature of 3 MeV. This corresponds to 
a very small zero sound oscillation. The initial conditions are varied in a very simple way, 
namely by changing the density p by an amount Ap = 10~ 4 po- No appreciable change of 
the results are obtained with different prescriptions of modifying the initial conditions, as 
discussed later. We checked that in all the calculations here presented the total energy 
was conserved up to 1% and that the particle number remains constant during the whole 
evolution. For a density equal to the saturation density p , fig. 1 (a), the density profiles 
of nuclear matter corresponding to different initial conditions keep close to each other for 
the whole considered time interval. This is typical of the "regular" region of a dynamical 
system, and indicates the stability of the dynamics with respect to small perturbations. On 
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the contrary, if the density p is chosen well inside the spinodal region, in this case p = ~po> 
the density profiles, after a short time, differ completely between each other, even if initially 
they are close within a part over 10 4 . This is the typical behaviour of a dynamical system in 
a chaotic region, where very small initial perturbations are rapidly amplified. The relevant 
feature of fig. 1 (b) is that the two density profiles do not differ simply in their amplitudes, 
as it would occur if the linear response regime would be valid, but mainly in their shape, 
which can occur only in a regime of strong non-linearity. We have checked that the results 
are stable with respect to the increase of the precision in the numerical integration. In fact, 
a decrease of the time integration step does not affect the results, and therefore the diver- 
gence in shape of the trajectories is an intrinsic dynamical property of the nuclear mean 
field in this region. The chaotic dynamics is intimately related with the strong non-linearity 
of the evolution, which implies a strong coupling among the different modes. This can be 
checked by Fourier transforming the density profiles at different times. We found that the 
different wave numbers increase their mixing more and more as the evolution proceeds. At 
the latest stage a wide range of wave numbers result to be excited, despite the fact that 
only one wave number k is initially populated. This point is illustrated in fig. 2 for a single 
trajectory corresponding to the full curve of fig.l (b). The final wave number spectral den- 
sity has therefore little resemblance with the initial one, and this feature was found to be 
independent from the particular initial k. In a linear regime, on the contrary, the system 
would be equivalent to a set of independent harmonic oscillators, one for each wave number 
k, with an imaginary frequency, corresponding to the unstable region. Such a system is 
integrable, and the corresponding flow in phase space is regular. It is important to note 
that the mixing takes place in an interval of time of about 20 fm/c, which is short with 
respect to a typical collision time or expansion time, characteristic of heavy ion collisions at 
intermediate energy. The calculations here presented are in two dimensions, but the char- 
acteristic times of chaotic dynamics are expected to be little dependent on dimensionality. 
A behaviour similar to that shown in fig.l (b) was found for a single trajectory in the case 
of ref. H where a dissipative term of Langevin type has been added to the Vlasov equation. 
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In principle a sensitive dependence on the initial conditions should be present also in that 
case. Calculations in that direction are in progress. To perform a more quantitative analysis, 
we follow the standard procedure in the theory of dynamical systems of characterizing the 
divergence among the phase space trajectories by a set of Lyapunov exponents |J. Let d(t) 
be the distance between two trajectories, along an unstable direction in phase space, at a 
given time t, and define 

x(t) = , (2) 

where do = d(0). Then the corresponding Lyapunov exponent is obtained by the limit 

Aqo = lim lim Aft) (3) 

provided the limit exists. The very definition of eq. (2) implies that asymptotically for large 
time and small enough perturbation the divergence is essentially exponential, namely that 
the equation of motion for the distance between the two trajectories is a linear one. This 
by no means implies that each trajectory follows linear equations of motion. In numerical 
applications one has to select a series of small values of do, a set of increasing times t and 
check that the corresponding values of \(t) stabilize around a definite value Aoo. In the 
present case the actual number of degrees of freedom and therefore the number of Lyapunov 
exponents is just equal to the number of cells in which the phase space is divided. In order 
to calculate an average value of the Lyapunov exponent, it is convenient to choose a distance 
which is connected with a global description of the system. We have therefore chosen as 
distance between two trajectories the difference in norm between their density profiles 

d(t) = E \P?\t)-p?\t)\/N e , (4) 

i 

where the index i runs over the N c cells in ordinary space, and p[ , p\ are the densities in 
the cell i for the trajectories 1 and 2 respectively. This definition neglects possible differences 
which can occur in momentum space between the two trajectories. The definition of eq. (4) 
is sufficient for the present analysis. It should include the contribution of all the unstable 



modes, in particular the ones which dominate for large t. In fig. 3 (a) it is displayed the 
time evolution of X(t) for different values of p/po, taking an initial do = 10 _4 p. It was 
checked that for do values of this size, X(t) stabilizes and is independent of do- One can 
distinguish two well defined regimes. For densities p well inside the spinodal region, the 
time dependence of X(t) is quite weak, resulting in the almost straight lines appearing at 
the top of fig.3 (a). This behaviour unambiguously indicates that the Lyapunov exponent 
in this case is well defined and that the system evolves in a genuine chaotic way. As the 
density increases and approaches the upper limit of the spinodal region, the general trend 
starts to change. In this case X(t) becomes a decreasing function of t which approaches 
rapidly very small values, indicating that the trajectories actually do not diverge between 
each other. The general behaviour appearing in fig.3 (a) is typical of dynamical systems 
possessing a regular and a chaotic region ||. The actual value of X(t) can be considered a 
measure of the non-linearity and chaoticity of the system ||. It has to be also stressed that, 
when the dynamics is non-linear, the Lyapunov exponent has not to be confused with 
the rate of increase of the density fluctuations. In the present case we found that, for a given 
trajectory, the density fluctuations increase with a characteristic time || substantially longer 
than the divergence characteristic time 1/Aqo [0, by a factor ranging from 2 to 3, according 
to the density and the initial conditions. In fig.3 (b) we show the Lyapunov exponent Aqo 
calculated at different densities. For density values close to the spinodal boundary X(t) 
shows some oscillations even for the largest times considered. However in this region the 
values are very small and the general trend of fig.3 (b) does not depend on this uncertainty. 
The value of X^ has a flat maximum around the density p = 0.35po- The curve reported in 
fig.3 (b) is essentially independent from k, which is in sharp contrast with the linear regime 
behaviour, for which the imaginary frequency is mainly a linear function of k. It can be 
instructive to compare the characteristic times tmf — fr/Xoo, which defines the time scale of 
the divergence between mean field trajectories, with the single particle characteristic time 
T sp — h/Ep, being Ep the Fermi energy at the given density. This is done in table I for a set 
of densities. One can see that for densities p < 0.4po the divergence time is smaller than the 
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single particle time, and therefore in this region the notion itself of mean field ceases to have 
any validity. In this region it appears more appropriate to speak of a fully many particle 
chaotic behaviour |K|]. It has to be stressed that, despite the density is small, the mean 
distance between particles is only slightly larger than at normal density, and therefore the 
particles can still strongly interact. In central collisions between heavy ions at intermediate 
energy nuclear matter is expected to be compressed at the initial stage of the reaction, until 
a highly excited composite nuclear system is formed. Numerous calculations |§, based on 
kinetic equations as well as on molecular dynamic models Rill , show that the momentum 
distribution is close to spherical symmetry and not too far from thermodynamic equilibrium. 
In this first stage of the reaction the dynamics is insensitive to small modifications of the 
initial conditions. Once the maximum compression is reached, nuclear matter starts to 
expand and can merge into the spinodal region. At this point the dynamics should change 
dramatically. According to the results presented in this paper, chaotic dynamics sets in, and 
the corresponding large density fluctuations dominate the cluster formations. The degree of 
chaoticity and fluctuations depends of course on the details of the dynamics which brings 
the reaction inside the spinodal region. The ratio between the expansion characteristic time 
and tmf could indicate whether chaos has the time to develop. The multifragmentation 
regime should appear in the energy interval where the chaoticity and the corresponding 
fluctuations are strong enough to produce several IMF. A confirmation of this scenario 
could be obtained by means of numerical simulations of heavy ion reactions. Unfortunately, 
at the moment, the extension of the calculations presented in this paper to the case of 
collisions between heavy ions appears problematic even in two dimensions. Thus one should 
find the signature of chaotic dynamics using more phenomenological approaches. The study 
of fragment size fluctuations and their statistical properties ||12|| , in particular along the 
lines of the intermittency analysis |13| , can be one of the viable possibilities. However, the 
problem is still completely open and demands a careful analysis, both theoretically and 
experimentally. 
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TABLES 

TABLE I. Comparison between the characteristic times tmf = 1/A<x> for the mean field and 
T sp = H/Ef for the single particle at different values of the nuclear density. See text. 



P/Po 








Jm/c 


Jm/c 


0.7 


7.76 


1538. 


0.6 


9.06 


28.57 


0.5 


10.87 


13.34 


0.4 


13.59 


9.92 


0.3 


18.12 


9.06 


0.2 


27.17 


10.00 
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FIGURES 

FIG. 1. The evolution for two close trajectories (dotted and full curve) for the nuclear matter 
density profile is illustrated as a function of time. Two different cases are shown for values of the 
initial density, (a) p/po=l and (b) p/po=0.5. While in case (a) the infinitesimal initial difference 
keeps constant in time and it is hardly visible, in case (b) it grows exponentially. Note the different 
scale used in fig. 1(a) and 1(b). 

FIG. 2. Power spectrum in momentum space of the density profile shown in fig. 1(b) as full 
curve. 

FIG. 3. (a) The behaviour of X(t) is plotted vs. time for different values of the initial density. 
The constant value obtained for p/po=OA, 0.5, 0.6 allows to define the Lyapunov exponent Aqo 
and to classify the chaoticity of the mean field at T=3 MeV. At variance for p/po=0.7 X(t) is 
a decreasing function of t and the mean field has a regular linear evolution, (b) The Lyapunov 
exponent Aqo is displayed for different values of the initial density. See text. 
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